W metodzie implicite zwanej schematem alfa, służącej do rozwiązywania problemów niestacjonarnych, nasz operator różniczkowy opisujący "fizykę" symulowanego zjawiska, stosujemy do kombinacji liniowej wartości stanu poprzedniego i następnego:
\( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}(\alpha u_t+(1-\alpha)u_{t+1})) = {f_{t+1}} \).
Jeśli operator różniczkowy jest liniowy, możemy wtedy rozdzialić go na dwa człony
\( u_{t+1} -dt \alpha \mathcal{L}(u_{t+1}) = u_t+dt{f_{t+1}}+ dt(1-\alpha)\mathcal{L}(u_t) \)
zostawiając po lewej stronie naszą "niewiadomą" \( u_{t+1} \).
Zauważmy że dla wartości alfa=1 schemat ten jest schematem implicte Eulera (tak zwany backward Euler), dla wartości alfa=0 jest on schematem explicite (tak zwany schemat Eulera) a dla wartości alfa=1/2 jest on schematem Cranka-Nicolson (który też jest schemate implicite).
Schemat Eulera został po raz pierwszy opisany w 1768 roku w podręczniku Leonharda Eulera. Schemat Cranka-Nicolson został wynaleziony w 1947 roku przez Phyllis Nicolson, brytyjską matematyczkę i Johna Cranka, amerykańskiego matematyka i fizyka [1], [2] .
Następnie stosujemy solwer izogeometrycznej metody elementów skończonych.
Żeby opracować solwer wykorzystujący metodę implicite w izogeometrycznej metodzie elementów skończonych, musimy przetransformować sformułowanie silne w sformułowanie słabe.
Przemnażamy więc nasze równanie przez funkcje testujące
\( (u_{t+1},w) -dt \alpha(\mathcal{L}(u_{t+1}),w) = (u_t,w)+dt({f_{t+1}},w)+ dt(1-\alpha)(\mathcal{L}(u_t),w) \)
Do aproksymacji stanu naszego systemu w danej chwili czasowej użyjemy kombinacji liniowej funkcji B-spline. W tym celu wybieramy bazę dwuwymiarowych funkcji B-spline, określając wektory węzłów w kierunku osi układu współrzędnych. Dla ustalenia uwagi możemy wybrać dwuwymiarową bazę funkcji B-spline drugiego stopnia
\( \{B^x_{i,2}(x)B^y_{j,2}(y)\}_{i=1,...,N_x;j=1,...,N_y} u_{i,j } \)
Będą one stosowane do aproksymacja symulowanego pola skalarnego aktualnej chwili czasowej
\( u(x,y;t+1)\approx \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } B^x_i(x)B^y_j(y) \)
Podobnie do testowania użyjemy funkcji bazowych B-spline:
\( w(x,y) \in \{ B^x_k(x)B^y_l(y) \}_{k=1,...,N_x;l=1,...,N_y} w^{k,l } \)
Zakładając że operator różniczkowy opisujący "fizykę" jest liniowy, nasze równanie wygląda teraz następująco:
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } (B^x_i(x)B^y_j(y)-dt \alpha \mathcal{L}(B^x_i(x)B^y_j(y)),B^x_k(x)B^y_l(y))=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \). Na przykład dla problemu transportu ciepła mamy
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ (B^x_i(x)B^y_j(y)-dt \alpha \left(\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2 }\right),B^x_k(x)B^y_l(y))=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \).
Dzięki sformułowaniu słabemu możemy scałkować przez części
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ \left( (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) - dt \alpha(\frac{\partial (B^x_i(x)B^y_j(y))}{\partial x},\frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}) -dt* (\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y},\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y }) \right)=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \).
Dzięki strukturze produktu tensorowego funkcji B-spline oraz dzięki temu że pochodna w kierunku y z funkcji B-spline w kierunku osi x daje 0 (ponieważ funkcje te są stałe w kierunku osi y) i podobnie dla pochodnej w kierunku y z funkcji B-spline w kierunku osi x mamy
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ \left( (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) - dt \alpha (\frac{\partial B^x_i(x)}{\partial x}B^y_j(y),\frac{\partial B^x_k(x)}{\partial x}B^y_l(y)) -dt* (B^x_i(x)\frac{\partial B^y_j(y)}{\partial y},B^x_k(x)\frac{\partial B^y_l(y)}{\partial y }) \right)=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \).
Zauważmy że nasz układ równań do rozwiązania to
\( \left(M_x \otimes M_y-dt \alpha S_x \otimes M_y -dt* M_x \otimes S_y\right) u^{t+1} = F(t) \)
gdzie
\( \{ M_x \otimes M_y \}_{i,j,k,l}= \int B^x_i(x)B^x_k(x)B^y_j(y)B^y_l(y) = \int B^x_i(x)B^y_j(y)B^x_k(x)B^y_l(y) = {\bf M}_{i,j,k,l } \) to macierz masowa będąca iloczynem Kroneckera dwóch jednowymiarowych macierzy masowych,
\( \{ S_x \otimes M_y\}_{i,j,k,l} = \int \frac{\partial B^x_i(x)}{\partial x}\frac{\partial B^x_k(x)}{\partial x}B^y_j(y)B^y_l(y) = \int \frac{\partial B^x_i(x)}{\partial x}B^y_j(y)\frac{\partial B^x_k(x)}{\partial x }B^y_l(y) \) to iloczyn Kronecker jednowymiarowej macierzy sztywności i jednowymiarowej macierzy masowej, oraz
\( \{ M_x \otimes S_y \}_{i,j,k,l} = \int B^x_i(x) B^x_k(x) \frac{\partial B^y_j(y)}{\partial y} \frac{\partial B^y_l(y)}{\partial y } = \int B^x_i(x)\frac{\partial B^y_j(y)}{\partial y}B^x_k(x)\frac{\partial B^y_l(y)}{\partial y } \) to iloczyn Kroneckera jednowymiarowej macierzy masowej i jednowymiarowej macierzy sztywności.
Każdą z tych macierzy potrafimy sfaktoryzować w czasie liniowym stosując algorytm solwera zmienno-kierunkowego. Jednakże faktoryzacja ich sumy w czasie liniowym nie jest już możliwa. Jest to możliwe dopiero gdy wprowadzimy schemat kroków czasowych który umożliwia rozdzielenie macierzy na podmacierze w podkrokach czasowych tak aby koszt faktoryzacji pozostawał liniowy.
Schemat Peaceman'a-Reachford'a pozwala na rozdzielenie kroku czasowego na dwa podkroki
\( \left(M_x \otimes M_y-dt \alpha S_x \otimes M_y\right) u^{t+1/2} = F(t+1/2)+\left(dt (1-\alpha) M_x \otimes S_y\right) u^{t} \)
\( \left(M_x \otimes M_y-dt \alpha M_x \otimes S_y\right) u^{t+1} = F(t+1/2)+\left(dt (1-\alpha) S_x \otimes M_y\right) u^{t+1/2 } \)
gdzie możemy scalić macierze lewej strony w jedną macierz o strukturze produktu Kroneckera, która może zostać sfaktoryzowana w czasie liniowym za pomocą algorytmu solwera zmienno-kierunkowego
\( \left(M_x-dt \alpha S_x \right)\otimes M_y u^{t+1/2} = F(t+1/2)+\left(dt (1-\alpha) M_x \otimes S_y\right) u^{t} \)
\( M_x \otimes \left( M_y-dt \alpha S_y\right) u^{t+1} = F(t+1/2)+\left(dt (1-\alpha) S_x \otimes M_y\right) u^{t+1/2 } \)
Przyglądnijmy się na koniec jak wygląda operator prawej strony
Funkcja prawej strony reprezentuje tutaj zmiany wywołane przez siłe działającą na system w czasie kroku czasowego, plus dodatkowo zmiany wywołane przez fizykę zjawiska w poprzednim kroku czasowym
\( (u_{t+1},w) -dt (1-\alpha)(\mathcal{L}(u_{t+1}),w) = (u_t,B^x_k(x)B^y_l(y))+dt({f_{t+1}},B^x_k(x)B^y_l(y))+ dt(1-\alpha)(\mathcal{L}(u_t),B^x_k(x)B^y_l(y)) \)
W szczególności nasza prawa strona nie jest jednak bitmapą, ale projekcją sumy trzech elementów:

1. Stanu naszego układu w poprzedniej chwili czasowej \( (u_t,w)= \sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) \) (również przemnożonego przez funkcje testującą i scałkowanego po obszarze). Oczywiście do przedstawiania stanu naszego układu w poprzednim kroku czasowym również wykorzystujemy kombinacje liniową funkcji bazowych B-spline \( u(x,y;t)=u_t=\sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } B^x_i(x)B^y_j(y) \)
2. Zmian wywołanych przez "fizykę" symulowanego zjawiska, w czasie kroku czasowego. Zmiany te policzone są poprzez zastosowanie operatora różniczkowego reprezentującego modelowane zjawisko do stanu układu w chwili poprzedniej. Stan układu oczywiście reprezentowany jest przez kombinacje liniową funkcji B-spline. Tak więc operator różniczkowy stosowany jest do funkcji B-spline \( (dt (1-\alpha) \mathcal{L}(u_t),w)=\sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } ({\cal L} (B^x_i(x)B^y_j(y)),B^x_k(x)B^y_l(y)) \). Na przykład dla problemu transportu ciepła mamy \( \sum_{i=1,...,N_x;j=1,...,N_y} dt (1-\alpha) u^{t}_{i,j } (\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2},B^x_k(x)B^y_l(y)) \). Dzięki sformułowaniu słabemu możemy scałkować przez części \( \sum_{i=1,...,N_x;j=1,...,N_y} dt (1-\alpha) u^{t}_{i,j } \left( (\frac{\partial (B^x_i(x)B^y_j(y))}{\partial x}B^y_j(y),\frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y},B^x_k(x)\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y})\right) \) co dzięki strukturze produktu tensorowego funkcji B-spline, oraz dzięki temu że pochodna w kierunku y z funkcji B-spline zmiennej x daje 0, oraz pochodna w kierunku x z funkcji B-spline zmiennej y daje, dostajemy więc \( \sum_{i=1,...,N_x;j=1,...,N_y} dt (1-\alpha) u^{t}_{i,j } \left( (\frac{\partial B^x_i(x)}{\partial x}B^y_j(y),\frac{\partial B^x_k(x)}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial B^y_j(y)}{\partial y},B^x_k(x)\frac{\partial B^y_l(y)}{\partial y})\right) \). Z punktu widzenia "macierzy" człon ten daje nam \( dt (1-\alpha) (S_x \otimes M_x +M_x\otimes S_y) u_t \)
3. Zmiany wywołane przez siłę działającą na system w czasie kroku czasowego \( (f_{t+1},w)= (f_{t+1}(x,y),B^x_k(x)B^y_l(y)) \).